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Abstract 

This paper describes a finite volume computational thermo-fluid dynamics method to 
solve for Navier-Stokes equations in conjunction with energy equation and 
thermodynamic equation of state in an unstructured coordinate system. The system of 
equations have been solved by a simultaneous Newton-Raphson method and compared 
with several benchmark solutions. Excellent agreements have been obtained in each case 
and the method has been found to be significantly faster than conventional 
Computational Fluid Dynamic(CFD) methods and therefore has the potential for 
implementation in Multi-Disciplinary analysis and design optimization in fluid and 
thermal systems. The paper also describes an algorithm of design optimization based on 
Newton-Raphson method which has been recently tested in a turbomachinery application. 

Nomenclature 


Symbol 

Description 

A 

Area (in^) 

c L 

Flow Coefficient 

d,D 

Diameter (in) 

g 

Gravitational Acceleration (ft/ sec^) 

£c 

Conversion Constant (“ 32.174 lb-ft/lbf-sec^) 

Kf 

Flow Resistance Coefficient (/lbf-sec2/(lb-ft)2 ) 

Krot 

Non-dimensional Rotating Flow Resistance Coefficient 

L 

Length (in) 

m 

Mass Flow Rate (Ib/sec) 

P 

Pressure (lbf/ in^) 

Re 

Reynolds Number (Re = puD/p) 

r 

Radius (in) 

S 

Momentum Source (lbf) 

T 

Temperature (° F) 

u, V 

Velocity (ft/sec) 

X 

Coordinate Direction (in). Quality 

y 

Coordinate Direction 

Greek 

P 

Density (lb/ft^) 

e 

Angle Between Branch Flow Velocity Vector and Gravity Vector (deg), 
Angle Between Neighboring Branches for Computing Shear (deg) 

CO 

Angular Velocity (rad/sec) 

£ 

Absolute Roughness (in) 

e/D 

Relative Roughness 

P 

Viscosity ( Ib/ft-sec) 
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1.0 Introduction 


Computational fluid dynamic methods have been extensively used to analyze a new 
design or to investigate operational problems in Aerospace, Automobile, Nuclear and 
Process Industries. Finite volume method [1] has been the most popular among the 
available method to solve for multi-dimensional Navier-Stokes equations. Several 
commercial CFD codes [2-6] are based on finite volume method. Most CFD models 
require extensive computer resources to perform a single analysis. For example a typical 
flow calculation in a compressor blade requires several hours in a supercomputer and 
typical turnaround time is one week. Although there was a need. Multi-disciplinary 
analysis and design optimization for fluid and thermal systems have not been found very 
practical in industries. 

The SIMPLE (Semi-Implicit Pressure Linked Equation) algorithm [7] has been widely 
used in all finite volume based computer codes mentioned earlier. Although there are 
several improvements to the original algorithm, SIMPLE is still based on successive 
substitution method. Momentum and continuity equations are solved separately. A 
pressure correction equation is derived from continuity equation using a linearized 
momentum equation. Under-relaxation is often necessary to ensure convergence. It 
slows down computational speed significantly. On the other hand simultaneous solution 
of all governing equations eliminates the need for under-relaxation and convergence is 
obtained much faster than any method using a successive substitution technique. 
Moreover, most commercially available computer codes use a structured coordinate 
system to define the flow configuration. In industrial applications majority of the flow 
situations can not be easily modeled using a structured coordinate system. Therefore a 
need exists for an efficient numerical algorithm for solving Navier-Stokes equation in an 
unstructured co-ordinate system which can be utilized for multi-disciplinary design 
optimization and analysis. The purpose of this paper is to describe the extension of the 
Generalized Fluid System Simulation Program (GFSSP) [8], a network flow analysis 
code to account for multidimensional flow effects. GFSSP discretizes an irregular flow 
domain in nodes and branches. The scalar transport equations such as mass and energy 
conservation equations are solved at the nodes where as the vector properties such as 
momentum are solved in branches. GFSSP has been extensively applied in propulsive 
flows [9- 1 0] and it demonstrated a satisfactory numerical stability in wide range of flow 
situations. In this paper an attempt has been made to extend GFSSP’s solver to handle 
multi-dimensional flow and demonstrate its use for a design problem. 

2.0 Unstructured Finite Volume Grid 

The unstructured finite volume grid network for GFSSP is shown in Figure 1, which 
shows connectivity of five nodes with four branches. In this figure node-i is connected 
with four neighboring nodes (J = 1 to 4) . In structured coordinate systems the number of 
neighboring nodes are restricted to 2, 4 and 6 for one, two and three dimensional systems 
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respectively. On the other hand for an unstructured system, there is no such restriction on 
the number of neighboring nodes. The index k represents fluid species. 
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Figure 1: Schematic of Nodes and Branches for an Unstructured Finite Volume Grid 

3.0 Conservation Equations 

3.1 Mass Conservation Equation 

The mass conservation equation for the ith node can be represented by: 


n 


j = l 


m.. =0 
ij 


Equation 1 implies that the net mass flow from a given 
words, the total mass flow rate into a node is equal to 
node. 


0 ) 

node must equate to zero. In other 
the total mass flow rate out of the 


3.2 Momentum Conservation Equation 

The one dimensional form of momentum equation for every branch takes the following 
form: 


7 i ( u r u uH p r p j) A 


'c 

Inertia 


Pressure 
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Gravity 
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Friction 


Centrifugal Source 


Equation 2 represents the balance of fluid forces acting on a given branch. Inertia, 
pressure, gravity, friction and centrifugal forces are considered in the conservation 
equation. In addition to five forces, a source term S has been provided in the equation to 
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input pump characteristics or to input power to pump in a given branch. If a pump is 
located in a given branch, all other forces except pressure are set to zero. The source term 
S is set to zero in all other cases. Figure 2 shows the schematic of a branch control 
volume. 
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Figure 2: Schematic of a Branch Control Volume Showing the Gravity and Rotation 

Multi-dimensional conservation equations must account for the transport of mass, 
momentum and energy into and out of the control volume from all directions in space. 
Mass conservation equations (Equation 1) can account for such transport because each 
internal node can be connected with multiple neighboring nodes located in space in any 
arbitrary location (Figure 1). On the other hand, the momentum conservation equation 
(Equation 2) is one dimensional. Multi -dimensional momentum transport can be 
accounted for by incorporating two additional terms in the momentum equation. These 
terms include: a) momentum transport due to shear, and b) momentum transport due to 
the transverse component of velocity. 

These two terms can be identified in the two-dimensional steady state Navier-Stokes 
equation which can be expressed as: 


du 

du^ 

5(puu) 

d(pvu) 

dp 

( d 2 u d 2 u > j 

v U dx 

+ v — 
d yj 

' dx 

II 

4 - 

-^ +p8 ‘ + * 

Idx 2 + dy 2 J 


The first term on the left hand side of Equation 3 corresponds to the current inertia term 
described in Equation 2 . The second term on the left hand side of Equation 3 
corresponds to the transverse momentum exchange. The first term on the right hand side 
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corresponds to the pressure term in Equation 2. The second term on the right hand side of 
Equation 3 corresponds to the gravity term. The third term on the right hand side of the 
equation is negligible (based on an order of magnitude argument). The fourth term on the 
right hand side represents momentum transport due to shear. This section discusses the 
implementation of shear and transverse momentum transport into the momentum 
equation of GFSSP. 


3.3 Momentum Transport Due to Shear 


Begin by examining the shear term (fourth term) of the Navier-Stokes Equation in more 
detail. First, consider the shear as a force instead of a force per unit volume by 
multiplying the volume by the shear term. 


5 2 u 

* P 

dy 


Au , ,/ Au, , Au 

^■(Ax)(Ay)(Az) = Axiz = 


shear 


( 4 ) 


Figure 3 represents a set of nodes and branches for which shear forces are exchanged. 
Let branch 1 2 represent the branch for which the shear force is to be calculated. Branches 
N12 and S12 represent the parallel branches which will be used to calculate the shear 
force on branch 12. Let YS be the distance between branches 12 and S12, and let YN be 
the distance between branches 12 and N12. Let AS be the shearing area between 
branches 12 and SI 2, while AN is the shearing area between branches 12 and N 12. 

A differencing scheme that can account for non-orthogonality in node structure was used. 
Equation 5 represents the shear term for branch 12. The angle 0 represents the angle that 
adjacent branches make with respect to the referenced branch. 





Branch 12 shear 



U NI2 COS ®N12 ~ U 12 

YN 



+ 



~ U S12 COS0S12 

YS 



( 5 ) 


Equation 5 can be generalized to n-number of parallel branches at any position around 
branch 12 as shown in Equation 6. 
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Figure 3: Branch and Node Schematic for Shear Exchange 


du ~ u il2 cos9 il2 -u 12 

P dy Branch 1 2 shear ~ “ v ^ Yi 

Where the summation subscript i represents the i‘ h parallel branch to branch 12. 



( 6 ) 


Next, the shear interaction with a neighboring parallel wall should be addressed. Suppose 
that adjacent to branch 12 is a wall that is approximately parallel to the branch. The angle 
Owaii represents the angle between branch 12 and the wall. The wall has a velocity v so , id . 
The distance between the centerline of branch 12 and the wall is y^,, and the shear area is 
A waM . The expression for the shear effect of the wall on branch 12 is given in Equation 7. 


3u 

dy 


p- 


U wall COS ® wall ~ U 12 


Branch 12 -wall^^. 


y wall 


L wall 


(7) 


If there were multiple walls adjacent to branch 12, then Equation 7 could be generalized 
into Equation 8. 
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( 8 ) 


Finally, combining the adjacent branch and the wall shear equations into one generalized 
equation for the i ,h branch for which shear is to be calculated. Equation 9 represents the 
actual laminar shear formulation incorporated into GFSSP, where i is the current branch, 
npj is the number of parallel branches to branch i, and ns, is the number of parallel solid 
walls to branch i. 
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3.4 Transverse Momentum Transport 


The transverse momentum component of Equation 3 can be expressed in terms of a force 
per unit volume. 

* M ( A x)( A y)( A z) = (A x)(A z)(pvAu) = m trans Au (10) 

dy (Ay) 

Figure 4 represents a set of nodes and branches for which transverse momentum 
exchange will take place. Let the branch 12 represent the current branch which will 
receive transverse momentum from the surrounding branches. Branch S12 represents an 
adjacent parallel branch, while branches SI and S2 represent the adjacent normal 
branches. 



Figure 4: Branch and Node Schematic for Transverse Momentum Exchange 


Now, examine the formulation for calculating the transverse momentum term for branch 
12. First, calculate the average mass flow rate for the adjacent normal branches: 


m s =^(m sl +m S2 ) 


( 11 ) 


Examining Figure 4, a positive transverse mass flow rate is defined as flow into the nodes 
corresponding to the branch in question. Based on this definition of transverse mass flow 
rate, calculate the transverse momentum term: 


5pvu 

dy 


V 


12 


(" 


m s u )2 


m s u sl2 


) = rh s (u 12 - u sl2 ) 


( 12 ) 


The parallel branch (S 1 2) will contribute to the transverse momentum term of branch 1 2 
when rh,, > 0 (since a positive transverse flow rate begins at S12 and ends at 12), and 
should have a negligible contribution when mi 2 <0 (i.e. transverse flow rate begins at 12 
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and ends at SI 2). Equation 13, below, is the upwinding representation of the transverse 
momentum term for branch 12. 


w u l2 (max||- m s ,o||) - u sl2 (max||o,m s ||) (13) 

12 

Equation 13 can be generalized for m-parallel branches around branch i, each with an 
angle 0^ with respect to branch i, and ny corresponding transverse connecting branches, 
each transverse branch with an angle of 0 ljk with respect to branch i. Equation 14 
represents this generalized version of Equation 13. 


Spvu 

dy 


V 


dpvu 

dy 


V 



(uu 


cosGijM max p 


k=i L n «i 


mijk cosBijk 


(14) 


Equation 14 represents the actual transverse momentum formulation put into GFSSP, 
where i is the current branch for which transverse momentum is being calculated, m, is 
the number of parallel branches which will be used to calculate transverse momentum, 
and n,j is the number of connecting transverse branches between the current branch (i), 
and the j ,h parallel branch. 


4.0 Verification Results 

In order to verify proper implementation of the shear and transverse momentum 
components into GFSSP, two models were identified and developed. The two 
verification models are: two dimensional Poiseuille flow and two dimensional shear 
driven flow in a square cavity. The following sections describe the models and present 
the results. 


4.1 Poiseuille Flow Model 

Consider the flow between two fixed flat plates shown below in Figure 5. The flow is 
pressure driven and assumed to be fully developed. The length of the plates is 1000 
inches, they are 1 inch apart, the fluid density is 12 lb m /ft\ the viscosity is 1 lb m /(ft sec), 
and the pressures at the upstream and downstream ends are 20 psi and 10 psi, 
respectively. The analytical solution (Equation 15) for this situation can easily be derived 
for laminar flow, and is provided below for comparative purposes. Figure 6 shows an 
approximate velocity profile for this situation. 
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Figure 5: Poiseulle Flow Physical Situation 



Figure 6: Poiseulle Flow Velocity Profile 


u = 0.005(y-y 2 ) 


(18) 


A simple 3 node, 10 branch model (two sets of 5 parallel branches) was constructed to 
model the physical situation described above. The model is shown schematically in 
Figure 7. 


Top Wall Top Wall 



Bottom Wall Bottom Wall 


Figure 7: GFSSP Poiseulle Flow Model 

Each of the 10 branches (121-125, 231-235) are 500 inches in length. It is assumed that 
the shear area for each branch-branch & branch-wall was 500 in 2 and that the shear 
distance was 0.2 inches between adjacent branches and 0.1 inch between walls and their 
adjacent branches. The pressures at nodes 1 and 3 are 20 psi and 10 psi, respectively. A 
constant density option was used and the corresponding density and viscosity used were 
1 21b m /ft 3 and 1.0 lb m /(ftsec), respectively. Resistance option -02 was used in the initial 
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flow field calculation (for the first Newton-Raphson iteration, after which the shear will 
replace the friction factor calculation) for each of the branches. The bottom and top walls 
are fixed. 

Figure 8 below shows a comparison between the velocity profiles for the analytical 
solution and the GFSSP 3 node, 10 branch (5 parallel branch) model. As can be seen in 
Figure 8, the results of this crude GFSSP model compare vary favorably with the 
analytical solution. 



0 0.1 0.2 0.3 0.4 0.5 0.6 

Velocity (ft/sec) 


Figure 8: Poiseulle Flow Velocity Distribution 


4.2 Shear Driven Flow in a Square Cavity 

Consider a square cavity as shown below in Figure 9. The flow is induced by shear 
interaction at the top wall. The length of each wall is 12 inches. The density of the fluid 
is assumed constant at 1.00 lb m /ft 3 , and the viscosity of the fluid is assumed to be 1.00 
lb m /(ft sec). The bottom and side walls are fixed. The top wall is moving to the right at 
constant speed of 100 ft/sec. The corresponding Reynolds number for this situation is Re 
= 100 . 


4.2.1 Benchmark Numerical Solution 
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Due to the non-linearity of the governing differential equations, an analytical solution of 
this situation is not available. Instead of an analytical solution, a well known numerical 
solution by Odus Burggraf [9] was used as the benchmark solution. Burggraf used a 
51x51 grid in his model of the square cavity. 



Figure 9: Flow in a Shear Driven Square Cavity 


4.2.2 GFSSP Driven Cavity Model 

The GFSSP model of the driven cavity consists of 49 nodes (48 of which are internal) 
and 84 branches. For numerical stability, one node (Node 1) was assigned to be a 
boundary node with an arbitrary pressure of 100 psi. A unit depth (1 inch) was assumed 
for the required areas. The model is shown schematically in Figure 10. 

It was assumed that the shear area for each branch is 2.0 in 2 . The shear distance between 
adjacent branches is 1.71429 inches and the shear distance between walls and their 
adjacent branches is 0.85714 inches. For transverse momentum, only adjacent parallel 
branches were considered, and only connecting branches not associated with the 
boundary node were used. A constant density option was used with corresponding 
density and viscosity values of 1.0 lb m /ft 3 and 1.0 lb m /(ft sec), respectively. Resistance 
Option 2 was used in the initial flow field calculation (for the first Newton-Raphson 
iteration, after which the shear will replace the friction factor calculation) for all of the 
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branches. The bottom and side walls are fixed. The top walls are moving to the right at 
100 ft/sec. All parallel angles are 0°, and all transverse angles are 90°. 
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Figure 10: GFSSP Flow in a Shear Driven Square Cavity Model 
4.2.3 Results 


Figure 1 1 shows a comparison between the benchmark numerical solution and GFSSP 
7x7 node model velocity profiles along a vertical plane at the horizontal midpoint. As 
can be seen in Figure 11, the results of this crude GFSSP model compare very favorably 
with the benchmark numerical solution of Burggraf [11]. 
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Dimensionless Velocity 

Figure 1 1 : Shear Driven Square Cavity Centerline Velocity Distribution 


5.0 Design Mode 

The design of a system is an iterative process, and all iterative schemes follow a ‘guess 
and correct’ procedure. With an intelligent guess one can reduce the required number of 
iterations and therefore can reduce the cost of the design of such system. In the past, 
every engineering organization depended on a few experienced engineers who can make 
such intelligent guesses by virtue of their long years of experience in a particular 
discipline. Because of the changing nature of the economy and technological 
developments, availability of such precious human resources are in decline. The future 
design of engineering hardware will be based on a computer based iterative process 
where an intelligent guess will no longer contribute to the economic aspect of the design. 
All design decisions will be based on computer models. 

5.1 Overview of the Design Mode Methodology 

This section describes how GFSSP’s analysis capability can be utilized for designing a 
fluid system. Suppose we want to specify a number of design variables in the flow 
circuit. These variables could be a) pressure at a given node(p„), b) flow rate in a 
particular branch(m b ) and c) temperature at a given node(T n ). It is necessary to identify 
the parameters that can be adjusted to attain the desired variables. Let us suppose these 
parameters are a) area of a particular branch, A x , b) Pressure at a given boundary node, p x , 
and c) Heat load at a given node, Q x . It should be noted that the number of desired 
variables must be equal to the number of identified parameters that can be adjusted during 
the iterative scheme. 
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Following are the seven steps of Newton-Raphson method proposed for design mode 
calculation: 

The application of the Newton-Raphson Method involves the following 7 steps: 

1 . Develop the governing equations. 

The equations are expressed in the following form: 

F \(P n >A x ,P x ,Q x ) = p„-f(A x ,p x ,Q x ) = o 

F 2^ b ,A x ,P x ,Q x ) = m b -f{A x ,P x ,Q x ) = 0 (16) 

■ F 3 ( T n ’ A X’P X ’Q x ) = mb-f( A x>P X ’Q x ) = ° 


2. Guess a solution for the equations. 

Guess A' x ,p‘ x ,Q[ as an initial solution for the governing equations 

3. Calculate the residuals of each equation. 

When the guessed solutions are substituted into Equation 16, the right hand side of the 
equation is not zero. The non-zero value is the residual. 

Fx(p n ,A\,p\,Q‘ x ) = p„-f ( A\,P X ,Q' X ) = R\ 

Fi(m b ,A' x ,p' x ,Q\) = m b -f ( A\,p\,Q\) = R 2 (17) 

Fi(T n ,A\,P x ,Q' x ) = T n -f ( A' x ,p x ,Q]) = R i 


The intent of the solution scheme is to correct A\,p x ,Q x with a set of corrections 

r t r 

A x ,p x ,Q x such that ^1,^2 >-^3 are zero. 


4. Develop a set of correction equations for all variables. 
First construct the matrix of influence coefficients: 
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( 18 ) 


0F . 

dF, 

dF, 

dA x 

d P x 

dQ x 

dF 2 

dF 2 

dF 2 

dA x 

d P x 


dF 3 

dF, 

ff 3 

dA x 

d P> 

FQ x 


Then construct the set of simultaneous equations for corrections: 


dF x 

dF x 

dF x 

FA X 

d P , 

dQ x 

dF 2 

dF i 

dF 2 

FA X 

d P , 

dQ x 

OF, 

FF 3 

FF 3 

dA x 

d P x 




(19) 


I t I 

5. Solve for A x ,p x ,Q x by solving the simultaneous equations. 


6. Apply correction to each variable. 

7. Iterate until the corrections become very small. 

5.2 An Example of the Design Mode Methodology 

5.2.1 Description 

It is required to estimate the cross-sectional area of a pressure regulator to set the 
downstream pressure at 55 psia. The pressure regulator is located 25” downstream in a 
(1/4” OD, 0.18” ID) pipeline with a relative roughness of 0.00006. The pipeline is 
connected with a 300 psia, 70° F nitrogen supply line. Nitrogen is vented to the 
atmosphere at 14.7 psia and 70° F at 60” downstream of the pressure regulator. The 
diameters and roughness are identical for both pipes. A flow restriction with a loss 
coefficient (C L ) of 0.6 is used to model the pressure regulator. Figure 12 schematically 
shows the situation to be modeled. 
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Atmosphere 
p = 14.7 psia 


Figure 12: Schematic of the Pressure Regulator System to be Modeled 


5.2.2 GFSSP Model 


The GFSSP model used to predict the pressure downstream of the pressure regulator is 
shown in Figure 13. Node 1 represents the tank, node 2 represents a position just 
upstream of the pressure regulator, node 3 represents a position just downstream of the 
regulator, and node 4 represents the atmosphere. Branch 23 represents the regulator, 
whereas branches 12 and 34 represent the pipes up- and downstream of the regulator. An 
initial guess of the cross-sectional area of the pressure regulator is made and a predicted 
pressure downstream of the regulator is calculated, which is likely to be different from 
the desired pressure of 55 psia. The intent of this exercise is to adjust the area until the 
desired pressure is achieved. The Design Mode Methodology, described in Section 5.1, 
is used adjust the regulator area within a convergence criterion of 0.01 on the desired 
pressure. 


-CH>PK 1> 3 


Figure 13: GFSSP Model of the Pressure Regulator System 


5.2.3 Results 


An initial guess of the pressure regulator cross-sectional area was made of 0.002 inches. 
Using this area, GFSSP predicted the pressure at node 3 to be 27.2 psia (well below the 
desired pressure of 55 psia). The Design Mode methodology, described in section 5.1, 
was employed to predict a new area for branch 23, and then a new predicted pressure was 
calculated for node 3. This process was repeated until the convergence criteria was met 
(three iterations were required). Table 1 summarizes the cycle of guessed and predicted 
values. 
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Table 1 : Summary of Guessed and Predicted values for Design Mode Example 


Iteration 

A 23 Guessed (in 2 ) 

P 3 Predicted (psia) 

Converged? 


1 

0.002 

27.20 

N 

0.005658 

2 

0.005658 

54.71 

N 

0.0057 

3 

0.0057 

55.00 

Y 



This systematic iterative procedure to adjust certain variables to achieve a design goal has 
also been demonstrated [12] in a turbomachinery application. 

6.0 Conclusions 

A finite volume based network analysis method was extended to perform multi- 
dimensional flow calculations in an unstructured co-ordinate system. Two classical flow 
problems were solved and results compared very well with benchmark analytical and 
numerical solution. Simultaneous solution of mass and momentum conservation 
equations was found more economical than conventional successive substitution 
employed in conventional CFD methods. The simultaneous method of solving Navier- 
Stokes equation make it more feasible to apply for multi-disciplinary analysis and design 
optimization. A design algorithm based on Newton-Raphson method was described and 
a simple application of the design method was demonstrated by solving a problem of 
finding the required area of a pressure regulator in a compressible flow feed line. 
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